library(data.table)
library(rmarkdown)
library(AnnotationHub)
library(tidyverse)
library(tximport)
library(ggplot2)
library(DESeq2)
library(pheatmap)
library(gridExtra)
library(ggplotify)
library(ggrepel)
library(UpSetR)
AnnotationSpecies <- "Homo sapiens" # Assign your species
ah <- AnnotationHub(hub=getAnnotationHubOption("URL")) # Bring annotation DB
ahQuery <- query(ah, c("OrgDb", AnnotationSpecies)) # Filter annotation of interest
if (length(ahQuery) == 1) {
DBName <- names(ahQuery)
} else if (length(ahQuery) > 1) {
DBName <- names(ahQuery)[1]
} else {
print("You don't have a valid DB")
rmarkdown::render()
}
AnnoDb <- ah[[DBName]] # Store into an OrgDb object
# Explore your OrgDb object with following accessors:
# columns(AnnpDb)
# keytypes(AnnoDb)
# keys(AnnoDb, keytype=..)
# select(AnnoDb, keys=.., columns=.., keytype=...)
AnnoKey <- keys(AnnoDb, keytype="ENSEMBLTRANS")
# Note: Annotation has to be done with not genome but transcripts
AnnoDb <- select(AnnoDb,
AnnoKey,
keytype="ENSEMBLTRANS",
columns="SYMBOL")
# Check if your AnnoDb has been extracted and saved correctely
class(AnnoDb)
## [1] "data.frame"
head(AnnoDb)
## ENSEMBLTRANS SYMBOL
## 1 ENST00000543404 A2MP1
## 2 ENST00000566278 A2MP1
## 3 ENST00000545343 A2MP1
## 4 ENST00000544183 A2MP1
## 5 ENST00000286479 NAT2
## 6 ENST00000520116 NAT2
.sf files have been created from fastq data by salmon
# This code chunk needs to be written by yourself
# Define sample names
SampleNames <- c("Mock_72hpi_S1",
"Mock_72hpi_S2",
"Mock_72hpi_S3",
"SARS-CoV-2_72hpi_S7",
"SARS-CoV-2_72hpi_S8",
"SARS-CoV-2_72hpi_S9")
# Define group level
GroupLevel <- c("Mock", "COVID")
# Define contrast for DE analysis
Contrast <- c("Group", "COVID", "Mock")
# Define a vector for comparing TPM vs Counts effect
TvC <- c("TPM", "Counts")
# Define .sf file path
sf <- c(paste0(SampleNames,
".salmon_quant/quant.sf"))
# Define sample groups
group <- c(rep("Mock", 3), rep("COVID", 3))
# Create metadata
metadata <- data.frame(Sample=factor(SampleNames, levels=SampleNames),
Group=factor(group, levels=GroupLevel),
Path=sf)
rownames(metadata) <- SampleNames
# Explore the metadata
print(metadata)
## Sample Group
## Mock_72hpi_S1 Mock_72hpi_S1 Mock
## Mock_72hpi_S2 Mock_72hpi_S2 Mock
## Mock_72hpi_S3 Mock_72hpi_S3 Mock
## SARS-CoV-2_72hpi_S7 SARS-CoV-2_72hpi_S7 COVID
## SARS-CoV-2_72hpi_S8 SARS-CoV-2_72hpi_S8 COVID
## SARS-CoV-2_72hpi_S9 SARS-CoV-2_72hpi_S9 COVID
## Path
## Mock_72hpi_S1 Mock_72hpi_S1.salmon_quant/quant.sf
## Mock_72hpi_S2 Mock_72hpi_S2.salmon_quant/quant.sf
## Mock_72hpi_S3 Mock_72hpi_S3.salmon_quant/quant.sf
## SARS-CoV-2_72hpi_S7 SARS-CoV-2_72hpi_S7.salmon_quant/quant.sf
## SARS-CoV-2_72hpi_S8 SARS-CoV-2_72hpi_S8.salmon_quant/quant.sf
## SARS-CoV-2_72hpi_S9 SARS-CoV-2_72hpi_S9.salmon_quant/quant.sf
# Assign sample names to the input (.sf) file path
names(sf) <- SampleNames
# Run tximport
# tpm vs original counts
# input sf: a factor of all .sf files' path
txi_tpm <- tximport(sf,
type="salmon",
tx2gene=AnnoDb,
countsFromAbundance="lengthScaledTPM", # Extracts TPM
ignoreTxVersion=T)
txi_counts <- tximport(sf,
type="salmon",
tx2gene=AnnoDb,
ignoreTxVersion=T)
# tpm
head(txi_tpm$counts)
## Mock_72hpi_S1 Mock_72hpi_S2 Mock_72hpi_S3 SARS-CoV-2_72hpi_S7
## A2MP1 5.077648 5.609035 0.8096504 0.00000
## A3GALT2 0.000000 0.000000 0.0000000 0.00000
## A4GNT 0.000000 2.974963 1.0038198 0.00000
## AACSP1 14.891510 19.954031 9.2490944 19.26746
## AADACL2 0.000000 0.000000 0.0000000 0.00000
## AADACL4 1.013374 0.000000 2.0033163 0.00000
## SARS-CoV-2_72hpi_S8 SARS-CoV-2_72hpi_S9
## A2MP1 0.00000 7.1022218
## A3GALT2 0.00000 0.0000000
## A4GNT 0.00000 0.9959273
## AACSP1 17.46992 4.5607796
## AADACL2 0.00000 0.0000000
## AADACL4 1.98163 0.0000000
dim(txi_tpm$counts)
## [1] 12202 6
# counts
head(txi_counts$counts)
## Mock_72hpi_S1 Mock_72hpi_S2 Mock_72hpi_S3 SARS-CoV-2_72hpi_S7
## A2MP1 6.186 7 1 0
## A3GALT2 0.000 0 0 0
## A4GNT 0.000 3 1 0
## AACSP1 16.000 22 10 21
## AADACL2 0.000 0 0 0
## AADACL4 1.000 0 2 0
## SARS-CoV-2_72hpi_S8 SARS-CoV-2_72hpi_S9
## A2MP1 0 3
## A3GALT2 0 0
## A4GNT 0 1
## AACSP1 10 5
## AADACL2 0 0
## AADACL4 2 0
dim(txi_counts$counts)
## [1] 12202 6
# Set a function creating dds and vsd
dds_vsd_fn <- function(txi) {
# Create a DESeq object (so-calledd dds)
des <- DESeqDataSetFromTximport(txi,
colData=metadata,
design=~Group)
# Create a vsd object (so-called vsd)
ves <- vst(des, blind=T)
# Output them as a list
return(list(dds=des, vsd=ves))
}
TPM <- dds_vsd_fn(txi_tpm)
Counts <- dds_vsd_fn(txi_counts)
# Outputs
# dds: TPM/Counts[[1]] or TPM/Counts[['dds']]
# vsd: TPM/Counts[[2]] or TPM/Counts[['vsd']]
# TPM
TPM[[1]]
## class: DESeqDataSet
## dim: 12202 6
## metadata(1): version
## assays(1): counts
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(0):
## colnames(6): Mock_72hpi_S1 Mock_72hpi_S2 ... SARS-CoV-2_72hpi_S8
## SARS-CoV-2_72hpi_S9
## colData names(3): Sample Group Path
head(counts(TPM[[1]]))
## Mock_72hpi_S1 Mock_72hpi_S2 Mock_72hpi_S3 SARS-CoV-2_72hpi_S7
## A2MP1 5 6 1 0
## A3GALT2 0 0 0 0
## A4GNT 0 3 1 0
## AACSP1 15 20 9 19
## AADACL2 0 0 0 0
## AADACL4 1 0 2 0
## SARS-CoV-2_72hpi_S8 SARS-CoV-2_72hpi_S9
## A2MP1 0 7
## A3GALT2 0 0
## A4GNT 0 1
## AACSP1 17 5
## AADACL2 0 0
## AADACL4 2 0
# Counts
Counts[[1]]
## class: DESeqDataSet
## dim: 12202 6
## metadata(1): version
## assays(2): counts avgTxLength
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(0):
## colnames(6): Mock_72hpi_S1 Mock_72hpi_S2 ... SARS-CoV-2_72hpi_S8
## SARS-CoV-2_72hpi_S9
## colData names(3): Sample Group Path
head(counts(Counts[[1]]))
## Mock_72hpi_S1 Mock_72hpi_S2 Mock_72hpi_S3 SARS-CoV-2_72hpi_S7
## A2MP1 6 7 1 0
## A3GALT2 0 0 0 0
## A4GNT 0 3 1 0
## AACSP1 16 22 10 21
## AADACL2 0 0 0 0
## AADACL4 1 0 2 0
## SARS-CoV-2_72hpi_S8 SARS-CoV-2_72hpi_S9
## A2MP1 0 3
## A3GALT2 0 0
## A4GNT 0 1
## AACSP1 10 5
## AADACL2 0 0
## AADACL4 2 0
# TPM
TPM[[2]]
## class: DESeqTransform
## dim: 12202 6
## metadata(1): version
## assays(1): ''
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(4): baseMean baseVar allZero dispFit
## colnames(6): Mock_72hpi_S1 Mock_72hpi_S2 ... SARS-CoV-2_72hpi_S8
## SARS-CoV-2_72hpi_S9
## colData names(4): Sample Group Path sizeFactor
# Counts
Counts[[2]]
## class: DESeqTransform
## dim: 12202 6
## metadata(1): version
## assays(1): ''
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(4): baseMean baseVar allZero dispFit
## colnames(6): Mock_72hpi_S1 Mock_72hpi_S2 ... SARS-CoV-2_72hpi_S8
## SARS-CoV-2_72hpi_S9
## colData names(3): Sample Group Path
using ‘avgTxLength’ from assays(dds), correcting for library size gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates
# Set a function estimating size factors, dispersions, and perform wald test
DESeqPrep_fn <- function(List) {
List[[1]] <- estimateSizeFactors(List[[1]])
List[[1]] <- estimateDispersions(List[[1]])
List[[1]] <- nbinomWaldTest(List[[1]])
return(List)
}
# Update dds with the function
Counts <- DESeqPrep_fn(Counts)
TPM <- DESeqPrep_fn(TPM)
sizeFactors(Counts[[1]])
## NULL
sizeFactors(TPM[[1]])
## Mock_72hpi_S1 Mock_72hpi_S2 Mock_72hpi_S3 SARS-CoV-2_72hpi_S7
## 1.1075601 1.0513174 0.8852784 1.3664601
## SARS-CoV-2_72hpi_S8 SARS-CoV-2_72hpi_S9
## 0.7661170 1.0051344
# Size factors don't exist in the Counts dds!
# Normalization factors are calculated in the Counts dds instead!
assays(Counts[[1]])
## List of length 6
## names(6): counts avgTxLength normalizationFactors mu H cooks
assays(TPM[[1]])
## List of length 4
## names(4): counts mu H cooks
colData(Counts[[1]])
## DataFrame with 6 rows and 3 columns
## Sample Group Path
## <factor> <factor> <character>
## Mock_72hpi_S1 Mock_72hpi_S1 Mock Mock_72hpi_S1.salmon..
## Mock_72hpi_S2 Mock_72hpi_S2 Mock Mock_72hpi_S2.salmon..
## Mock_72hpi_S3 Mock_72hpi_S3 Mock Mock_72hpi_S3.salmon..
## SARS-CoV-2_72hpi_S7 SARS-CoV-2_72hpi_S7 COVID SARS-CoV-2_72hpi_S7...
## SARS-CoV-2_72hpi_S8 SARS-CoV-2_72hpi_S8 COVID SARS-CoV-2_72hpi_S8...
## SARS-CoV-2_72hpi_S9 SARS-CoV-2_72hpi_S9 COVID SARS-CoV-2_72hpi_S9...
colData(TPM[[1]])
## DataFrame with 6 rows and 4 columns
## Sample Group Path
## <factor> <factor> <character>
## Mock_72hpi_S1 Mock_72hpi_S1 Mock Mock_72hpi_S1.salmon..
## Mock_72hpi_S2 Mock_72hpi_S2 Mock Mock_72hpi_S2.salmon..
## Mock_72hpi_S3 Mock_72hpi_S3 Mock Mock_72hpi_S3.salmon..
## SARS-CoV-2_72hpi_S7 SARS-CoV-2_72hpi_S7 COVID SARS-CoV-2_72hpi_S7...
## SARS-CoV-2_72hpi_S8 SARS-CoV-2_72hpi_S8 COVID SARS-CoV-2_72hpi_S8...
## SARS-CoV-2_72hpi_S9 SARS-CoV-2_72hpi_S9 COVID SARS-CoV-2_72hpi_S9...
## sizeFactor
## <numeric>
## Mock_72hpi_S1 1.107560
## Mock_72hpi_S2 1.051317
## Mock_72hpi_S3 0.885278
## SARS-CoV-2_72hpi_S7 1.366460
## SARS-CoV-2_72hpi_S8 0.766117
## SARS-CoV-2_72hpi_S9 1.005134
# Extract and save the size factors in a data frame
sizeFactor <- as.data.frame(round(sizeFactors(TPM[[1]]), 3))
colnames(sizeFactor) <- 'Size_Factor'
sizeFactor <- sizeFactor %>%
rownames_to_column(var="Sample") %>%
inner_join(metadata[, 1:ncol(metadata)-1], by="Sample")
# Create a plot comparing the size factors by sample
ggplot(sizeFactor, aes(x=Sample,
y=Size_Factor,
fill=Group,
label=Size_Factor)) +
geom_bar(stat="identity", width=0.8) +
theme_bw() +
ggtitle("Size Factors in TPM-DESeq") +
geom_text(vjust=1.5) +
theme(axis.text.x=element_text(angle=45,
vjust=0.5)) + ylab("Size Factor") + geom_hline(yintercept=1, color="blue", linetype="dashed")
# Extract and normalization factors in a data frame
normf <- as.data.frame(normalizationFactors(Counts[[1]])) %>%
gather(Sample, Normalization_Factor) %>%
inner_join(metadata[, 1:2], by="Sample")
normf$Sample <- factor(normf$Sample, levels=SampleNames)
normf$Group <- factor(normf$Group, levels=GroupLevel)
# Create a scatter plot showing distribution of normalization factors
normFactors_plot <- ggplot(normf,
aes(x=Sample, y=Normalization_Factor)) +
geom_jitter(alpha=0.5, aes(color=Group)) +
# Add a boxplot to provide statistics in each sample
geom_boxplot(aes(x=Sample, y=Normalization_Factor),
outlier.shape=NA, alpha=0.5) +
theme_bw() +
ggtitle("Normalization Factors in Counts-DESeq") +
theme(axis.text.x=element_text(angle=45,
vjust=0.5)) +
ylab("Normalization Factor / Gene") +
# Add a dashed horizontal line to indicate where normalization factor=1
geom_hline(yintercept=1,
color="blue",
linetype="dashed")
# Print the normalization factor scatter plot
print(normFactors_plot)
# Print the same plot with larger y-magnification in order to observe the box plot
normFactors_plot +
ylim(0.5, 1.5)
# Assigne what to compare
GroupOfInterest <- Contrast[1]
# Set a function for a PCA plot
QCPCA_fn <- function(inputList, Title) {
plotPCA(inputList[[2]], # takes vsd
intgroup=GroupOfInterest) + theme_bw() + ggtitle(Title)
}
# Set heatmap annotation
ColOfInterest <- !colnames(metadata) %in% c("Sample", "Path")
HeatmapAnno <- as.data.frame(metadata[, ColOfInterest])
rownames(HeatmapAnno) <- SampleNames
colnames(HeatmapAnno) <- colnames(metadata)[ColOfInterest]
# Set a function for a correlation heatmap
QCcorrHeatmap_fn <- function(inputList, Title) {
# Extract transformed count matrix
mtx <- assay(inputList[[2]]) # takes vsd
# Calculate correlation and store in the matrix
mtx <- cor(mtx)
# Create a correlation heatmap
return(pheatmap(mtx,
annotation=HeatmapAnno,
main=paste("Sample Correlation Heatmap:",
Title)))
}
grid.arrange(QCPCA_fn(TPM, "QC PCA: TPM"),
QCPCA_fn(Counts, "QC PCA: Counts"),
nrow=2)
# TPM
QCcorrHeatmap_fn(TPM, "TPM")
# Counts
QCcorrHeatmap_fn(Counts, "Counts")
# Create a list for TPM and Counts dds
ddsList <- list(TPM=TPM[[1]], Counts=Counts[[1]])
for (x in TvC) {
# Run DESeq()
ddsList[[x]] <- DESeq(ddsList[[x]])
print(resultsNames(ddsList[[x]]))
}
## [1] "Intercept" "Group_COVID_vs_Mock"
## [1] "Intercept" "Group_COVID_vs_Mock"
# Plot dispersion
for (x in TvC) {
plotDispEsts(ddsList[[x]],
ylab="Dispersion",
xlab="Mean of Normalized Counts",
main=paste("Dispersion of", x, "Input"))
}
# Create an empty list for shrunken dds
shr_ddsList <- list(TPM=c(), Counts=c())
for (x in TvC) {
# shrink
shr_ddsList[[x]] <- lfcShrink(ddsList[[x]],
contrast=Contrast, # contrast
type="normal") # is paired with "normal" type
}
# Set FDR threshold
alpha=0.1
# FDR threshold vector
FDRv=c("< 0.1", "> 0.1")
# Set a function cleaning table
Sig_fn <- function(df, Input) {
df <- df %>%
rownames_to_column(var="Gene") %>%
mutate(FDR=ifelse(padj < 0.1 & !is.na(padj),
FDRv[1],
FDRv[2]),
Input=Input)
return(df)
}
# Initialize lists of result tables with (resList) or without (shr_resList) shrinkage
resList <- ddsList
shr_resList <- ddsList
for (x in TvC) {
# Extract results
resList[[x]] <- as.data.frame(results(ddsList[[x]],
contrast=Contrast,
alpha=alpha))
shr_resList[[x]] <- as.data.frame(shr_ddsList[[x]])
# clean the data frame
resList[[x]] <- Sig_fn(resList[[x]], x)
shr_resList[[x]] <- Sig_fn(shr_resList[[x]], x)
}
# No shrinkage summary
summary(resList)
## Length Class Mode
## TPM 9 data.frame list
## Counts 9 data.frame list
head(resList[['TPM']])
## Gene baseMean log2FoldChange lfcSE stat pvalue padj
## 1 A2MP1 3.0525639 -0.71117275 1.7535703 -0.4055570 0.6850681 NA
## 2 A3GALT2 0.0000000 NA NA NA NA NA
## 3 A4GNT 0.8296737 -1.91297083 3.0641144 -0.6243144 0.5324211 NA
## 4 AACSP1 13.9670247 -0.07686924 0.7514137 -0.1022995 0.9185190 0.9811568
## 5 AADACL2 0.0000000 NA NA NA NA NA
## 6 AADACL4 0.9621048 -0.34400388 2.9069797 -0.1183372 0.9058005 NA
## FDR Input
## 1 > 0.1 TPM
## 2 > 0.1 TPM
## 3 > 0.1 TPM
## 4 > 0.1 TPM
## 5 > 0.1 TPM
## 6 > 0.1 TPM
head(resList[['Counts']])
## Gene baseMean log2FoldChange lfcSE stat pvalue padj
## 1 A2MP1 2.8317018 -0.8410796 1.8117650 -0.4642322 0.6424814 NA
## 2 A3GALT2 0.0000000 NA NA NA NA NA
## 3 A4GNT 0.8395428 -1.9116000 3.0121489 -0.6346300 0.5256698 NA
## 4 AACSP1 13.9323130 -0.1237197 0.7552029 -0.1638231 0.8698704 0.9661114
## 5 AADACL2 0.0000000 NA NA NA NA NA
## 6 AADACL4 0.9757540 -0.3681081 2.8561649 -0.1288819 0.8974511 NA
## FDR Input
## 1 > 0.1 Counts
## 2 > 0.1 Counts
## 3 > 0.1 Counts
## 4 > 0.1 Counts
## 5 > 0.1 Counts
## 6 > 0.1 Counts
# Shrinkage summary
summary(shr_resList)
## Length Class Mode
## TPM 9 data.frame list
## Counts 9 data.frame list
head(shr_resList[['TPM']])
## Gene baseMean log2FoldChange lfcSE stat pvalue padj
## 1 A2MP1 3.0525639 -0.13588716 0.3372049 -0.4055570 0.6850681 NA
## 2 A3GALT2 0.0000000 NA NA NA NA NA
## 3 A4GNT 0.8296737 -0.12996851 0.2230905 -0.6243144 0.5324211 NA
## 4 AACSP1 13.9670247 -0.04328458 0.4236828 -0.1022995 0.9185190 0.9811568
## 5 AADACL2 0.0000000 NA NA NA NA NA
## 6 AADACL4 0.9621048 -0.02659239 0.2313921 -0.1183372 0.9058005 NA
## FDR Input
## 1 > 0.1 TPM
## 2 > 0.1 TPM
## 3 > 0.1 TPM
## 4 > 0.1 TPM
## 5 > 0.1 TPM
## 6 > 0.1 TPM
head(shr_resList[['Counts']])
## Gene baseMean log2FoldChange lfcSE stat pvalue padj
## 1 A2MP1 2.8317018 -0.14479543 0.3333802 -0.4642322 0.6424814 NA
## 2 A3GALT2 0.0000000 NA NA NA NA NA
## 3 A4GNT 0.8395428 -0.13464596 0.2270494 -0.6346300 0.5256698 NA
## 4 AACSP1 13.9323130 -0.06903748 0.4245107 -0.1638231 0.8698704 0.9661114
## 5 AADACL2 0.0000000 NA NA NA NA NA
## 6 AADACL4 0.9757540 -0.02945511 0.2355359 -0.1288819 0.8974511 NA
## FDR Input
## 1 > 0.1 Counts
## 2 > 0.1 Counts
## 3 > 0.1 Counts
## 4 > 0.1 Counts
## 5 > 0.1 Counts
## 6 > 0.1 Counts
# Set ylim: has to adjusted by users depending on data
yl <- c(-20, 20)
# Set min log2 fold change of interest
mLog <- c(-1, 1)
# Define a function creating an MA plot
MA_fn <- function(List, Shr) {
MAList <- ddsList
for (i in 1:2) {
MAplot <- ggplot(List[[i]],
aes(x=baseMean,
y=log2FoldChange,
color=FDR)) + geom_point() + scale_x_log10() + theme_bw() + scale_color_manual(values=c("blue", "grey")) + ggtitle(paste("MA plot:", names(List)[i], "Input with", Shr)) + ylim(yl[1], yl[2])+ geom_hline(yintercept=c(mLog[1], mLog[2]), linetype="dashed", color="red")
MAList[[i]] <- MAplot
}
return(MAList)
}
# Create MA plots with or without shrinkage and store in a list
MA <- MA_fn(resList, "No Shrinkage")
shr_MA <- MA_fn(shr_resList, "Shrinkage")
# TPM with or without shrinkage
grid.arrange(MA[[1]], shr_MA[[1]], nrow=1)
# TPM with or without shrinkage
grid.arrange(MA[[2]], shr_MA[[2]], nrow=1)
# Combining total data table
res <- rbind(shr_resList[['TPM']], shr_resList[['Counts']])
res$Input <- factor(res$Input, levels=TvC) # TvC=c("TPM", "Counts")
# Create a plot presenting distribution of FDR
ggplot(res,
aes(x=padj, color=Input)) +
geom_density(size=1, aes(y=..count..)) +
theme_bw() +
ggtitle("Distribution of False Discovery Rate (FDR)") +
xlab("Adjusted P-Value (FDR)") +
ylab("Count") +
geom_vline(xintercept=alpha,
color="black",
linetype="dashed",
size=1) +
scale_x_continuous(breaks=seq(0, 1, by=0.1))
# Set xlim for volcano plots
xlim=c(-6, 6) # has to be assined by users
# Set a basic volcano plot function
Volcano_fn <- function(df, Label=NULL) {
ggplot(df,
aes(x=log2FoldChange,
y= -log10(padj),
color=FDR,
label=Label)) +
geom_point() +
facet_grid(.~Input) +
theme_bw() +
scale_color_manual(values=c("blue", "grey")) +
ggtitle("Volcano Plot") +
ylab("-log10(padj)") +
theme(strip.text.x=element_text(size=12)) +
geom_vline(xintercept=c(mLog[1], mLog[2]),
color="red",
linetype="dashed",
size=1) +
xlim(xlim[1], xlim[2])
}
# Display the volcano plots by input
Volcano_fn(res)
# Assign log odds threshold
LogOddsCut=20
# Add a column indicating high log odds genes
res <- res %>%
mutate(Label=ifelse(-log10(padj) > LogOddsCut,
Gene,
""))
# Display the genes with volcano plots
Volcano_fn(res, Label=res$Label) + geom_text_repel(color="black")
ggplot(res[res$FDR == "< 0.1", ], # Subset rows with FDR < alpha
aes(x=log2FoldChange,
color=Input)) +
geom_density(size=1, aes(y=..count..)) +
theme_bw() +
ylab("Count") +
geom_vline(xintercept=c(mLog[1], mLog[2]),
color="black",
linetype="dashed", size=1) +
ggtitle("Distribution of Log2 Folds by Input Type") +
xlim(xlim[1], xlim[2])
# Initialize a list
lowfdrList <- ddsList # A list for normalized counts matrix with FDR below alpha
highfoldList <- ddsList # A list for normalized counts with log2foldchange over minLog
for (x in TvC) {
# Create filtering vectors with alpha and log2foldchange
BelowAlpha <- shr_resList[[x]]$FDR == FDRv[1]
overmLog <- abs(shr_resList[[x]]$log2FoldChange) > mLog[2] # mLog has been set to c(-1, 1) previously
# Extract transformed counts from vsd objects (TPM[['vsd']] and Counts[['vsd']])
if (x == "TPM") {
normCounts <- counts(TPM[['dds']], normalized=T)
} else {
normCounts <- counts(Counts[['dds']], normalized=T)
}
# Update the normalized count matrix with FDR below alpha
lowfdrList[[x]] <- normCounts[BelowAlpha, ]
highfoldList[[x]] <- normCounts[BelowAlpha & overmLog, ]
summary(lowfdrList[[x]])
summary(highfoldList[[x]])
}
# Initialize map lists
lowfdrMap <- ddsList
highfoldMap <- ddsList
# Set a function creating a heatmap
ProfileHeatmap_fn <- function(inputmatrix, Title1, Title2, Title3=NULL) {
as.ggplot(pheatmap(inputmatrix,
annotation=HeatmapAnno,
scale="row", # presents z-score instead of counts
show_rownames=F,
main=paste("Transcription Profiles with",
Title1,
"input and",
Title2,
alpha,
Title3)))
}
# Create and save heatmaps
for (x in TvC) {
lowfdrMap[[x]] <- ProfileHeatmap_fn(lowfdrList[[x]],
Title1=x,
Title2="FDR <")
highfoldMap[[x]] <- ProfileHeatmap_fn(highfoldList[[x]],
Title1=x,
Title2="FDR <",
Title3=paste("+ Absolte Log2 Fold Change >", mLog[2]))
}
# Count number of NA genes
type=c("Zero Counts", "Outliers", "Total NA Genes")
NAstat <- res %>%
group_by(Input) %>%
summarize(zero=sum(is.na(log2FoldChange)),
outlier=sum(is.na(pvalue) & is.na(padj))) %>%
mutate(total=sum(zero, outlier)) %>%
gather(Type, Number, -Input) %>%
mutate(Type=factor(case_when(Type == "zero" ~ type[1],
Type == "outlier" ~ type[2],
Type == "total" ~ type[3]),
levels=type))
# Plot number of NA genes
ggplot(NAstat, aes(x=Type, y=Number, group=Input, fill=Input, label=Number)) +
geom_bar(stat="identity", position="dodge") +
theme_bw() +
geom_text(position=position_dodge(width=1), vjust=1.5) +
ggtitle("Number of NA Genes") +
ylab("Number of Genes")
# Create a new list having DE table with FDR below alpha
lowfdr_resList <- shr_resList
for (x in TvC) {
lowfdr_resList[[x]] <- filter(shr_resList[[x]],
FDR == FDRv[1]) %>%
as.data.table()
}
# Initialize new lists in order to store rank-updated result DE tables
FDRrankList <- lowfdr_resList
lfcList <- lowfdr_resList
UPlfcList <- lowfdr_resList
DOWNlfcList <- lowfdr_resList
# Set a function creating a column for the rank
Ranking_fn <- function(x) {mutate(x, Rank=1:nrow(x))}
for (x in TvC) {
# Rearrange genes with FDR
FDRrankList[[x]] <- lowfdr_resList[[x]][order(padj),] %>%
Ranking_fn()
# Rearrange genew with absolute log2FoldChange
lfcList[[x]] <- lowfdr_resList[[x]][order(-abs(log2FoldChange)),] %>%
Ranking_fn()
# Rearrange genes with log2FoldChange (decreasing order)
UPlfcList[[x]] <- lowfdr_resList[[x]][order(-log2FoldChange),] %>%
Ranking_fn()
# Rearrange genes with log2FoldChange (increasing order)
DOWNlfcList[[x]] <- lowfdr_resList[[x]][order(log2FoldChange),] %>%
Ranking_fn()
}
# Explore the ranks
print(c(FDRrankList, lfcList, UPlfcList, DOWNlfcList))
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: CCN2 3782.15461 3.3138239 0.1852169 17.879269 1.710657e-71
## 2: F2R 2460.48267 1.9146387 0.1209837 15.821943 2.196087e-56
## 3: CDCP1 1500.91438 2.6204056 0.1822811 14.363152 8.812579e-47
## 4: STK38L 11846.66658 1.8911153 0.1338459 14.128446 2.536811e-45
## 5: TM4SF18 1895.68247 1.9727888 0.1415730 13.930193 4.152312e-44
## ---
## 519: MRPL51 1960.13310 0.3127969 0.1272269 2.458574 1.394902e-02
## 520: L3MBTL3 291.07491 -0.4897894 0.1994004 -2.456223 1.404062e-02
## 521: SDHAF3 2795.76859 -0.3085702 0.1257367 -2.454098 1.412385e-02
## 522: HEXIM1 337.54914 -0.4614688 0.1882160 -2.451750 1.421635e-02
## 523: TMSB15B 77.90622 0.6931845 0.2827217 2.449988 1.428610e-02
## padj FDR Input Rank
## 1: 6.249030e-68 < 0.1 TPM 1
## 2: 4.011153e-53 < 0.1 TPM 2
## 3: 1.073078e-43 < 0.1 TPM 3
## 4: 2.316742e-42 < 0.1 TPM 4
## 5: 3.033679e-41 < 0.1 TPM 5
## ---
## 519: 9.818068e-02 < 0.1 TPM 519
## 520: 9.863534e-02 < 0.1 TPM 520
## 521: 9.902957e-02 < 0.1 TPM 521
## 522: 9.948720e-02 < 0.1 TPM 522
## 523: 9.978415e-02 < 0.1 TPM 523
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: CCN2 3838.81312 3.3117030 0.1858269 17.809254 5.990320e-71
## 2: F2R 2498.11791 1.9128696 0.1211696 15.783102 4.067010e-56
## 3: CDCP1 1523.30246 2.6177630 0.1826892 14.316349 1.729741e-46
## 4: STK38L 12016.77780 1.8894025 0.1345047 14.046514 8.092088e-45
## 5: TM4SF18 1919.12630 1.9716534 0.1416724 13.912413 5.325183e-44
## ---
## 527: CUTALP 80.56840 -0.6867126 0.2798177 -2.451451 1.422816e-02
## 528: SDHAF3 2839.34437 -0.3094651 0.1263165 -2.449916 1.428894e-02
## 529: CREG1 273.30711 -0.4547095 0.1859437 -2.445053 1.448308e-02
## 530: P2RY13 74.69163 -0.6867880 0.2807198 -2.444416 1.450869e-02
## 531: RBM43 33.32318 -0.8843602 0.3610070 -2.444555 1.450309e-02
## padj FDR Input Rank
## 1: 2.188863e-67 < 0.1 Counts 1
## 2: 7.430427e-53 < 0.1 Counts 2
## 3: 2.106825e-43 < 0.1 Counts 3
## 4: 7.392122e-42 < 0.1 Counts 4
## 5: 3.891644e-41 < 0.1 Counts 5
## ---
## 527: 9.867185e-02 < 0.1 Counts 527
## 528: 9.888595e-02 < 0.1 Counts 528
## 529: 9.983948e-02 < 0.1 Counts 529
## 530: 9.983948e-02 < 0.1 Counts 530
## 531: 9.983948e-02 < 0.1 Counts 531
##
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: EDN1 141.58628 3.4835855 0.3262539 10.347843 4.280208e-25
## 2: CCN2 3782.15461 3.3138239 0.1852169 17.879269 1.710657e-71
## 3: GPR87 526.04952 3.2831219 0.2391306 13.655049 1.883633e-42
## 4: MUC12 31.29951 3.0505542 0.4135533 5.990845 2.087531e-09
## 5: SERPINE1 8891.97515 2.8509705 0.3263039 8.759388 1.963139e-18
## ---
## 519: PPIL4 3537.05463 0.3036901 0.1171944 2.591328 9.560620e-03
## 520: ND1 35591.46112 -0.2990627 0.1045701 -2.859924 4.237422e-03
## 521: EIF4EBP2 4471.66924 0.2985887 0.1171230 2.549363 1.079198e-02
## 522: SF3B6 4172.24628 0.2930706 0.1117374 2.622851 8.719732e-03
## 523: RPL36AP37 15.01076 0.1506778 0.1855742 5.300249 1.156447e-07
## padj FDR Input Rank
## 1: 1.563560e-22 < 0.1 TPM 1
## 2: 6.249030e-68 < 0.1 TPM 2
## 3: 1.146819e-39 < 0.1 TPM 3
## 4: 1.229960e-07 < 0.1 TPM 4
## 5: 3.585673e-16 < 0.1 TPM 5
## ---
## 519: 7.446683e-02 < 0.1 TPM 519
## 520: 3.989511e-02 < 0.1 TPM 520
## 521: 8.128477e-02 < 0.1 TPM 521
## 522: 6.985347e-02 < 0.1 TPM 522
## 523: 4.912212e-06 < 0.1 TPM 523
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: EDN1 143.61165 3.4934078 0.3252202 10.406217 2.322691e-25
## 2: CCN2 3838.81312 3.3117030 0.1858269 17.809254 5.990320e-71
## 3: GPR87 533.54706 3.2823931 0.2390153 13.656904 1.836253e-42
## 4: ALPK2 44.93538 3.1746617 0.4192791 5.808344 6.309395e-09
## 5: MUC12 19.79726 2.9517456 0.4237826 6.511310 7.449834e-11
## ---
## 527: PPIL4 3591.69222 0.3024157 0.1176344 2.570801 1.014635e-02
## 528: ND1 36140.65876 -0.3001565 0.1054619 -2.846112 4.425659e-03
## 529: EIF4EBP2 4540.89362 0.2974076 0.1179816 2.520798 1.170890e-02
## 530: SF3B6 4236.85234 0.2917483 0.1123824 2.596033 9.430688e-03
## 531: RPL36AP37 15.15450 0.1511406 0.1861197 5.298424 1.168063e-07
## padj FDR Input Rank
## 1: 8.487113e-23 < 0.1 Counts 1
## 2: 2.188863e-67 < 0.1 Counts 2
## 3: 1.118278e-39 < 0.1 Counts 3
## 4: 3.440975e-07 < 0.1 Counts 4
## 5: 5.917760e-09 < 0.1 Counts 5
## ---
## 527: 7.723704e-02 < 0.1 Counts 527
## 528: 4.094014e-02 < 0.1 Counts 528
## 529: 8.660795e-02 < 0.1 Counts 529
## 530: 7.316291e-02 < 0.1 Counts 530
## 531: 4.962910e-06 < 0.1 Counts 531
##
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: EDN1 141.58628 3.483585 0.3262539 10.347843 4.280208e-25
## 2: CCN2 3782.15461 3.313824 0.1852169 17.879269 1.710657e-71
## 3: GPR87 526.04952 3.283122 0.2391306 13.655049 1.883633e-42
## 4: MUC12 31.29951 3.050554 0.4135533 5.990845 2.087531e-09
## 5: SERPINE1 8891.97515 2.850970 0.3263039 8.759388 1.963139e-18
## ---
## 519: TRMT9B 145.07064 -1.716032 0.3049279 -5.610108 2.021999e-08
## 520: TBC1D3L 28.55182 -1.742442 0.3813485 -4.499991 6.795649e-06
## 521: CTXND2 168.84610 -1.817134 0.2652655 -6.827901 8.616600e-12
## 522: WFIKKN2 371.00721 -1.937066 0.2298615 -8.412662 4.008097e-17
## 523: MROH7-TTC4 19.77094 -2.336210 0.4243946 -5.001629 5.684777e-07
## padj FDR Input Rank
## 1: 1.563560e-22 < 0.1 TPM 1
## 2: 6.249030e-68 < 0.1 TPM 2
## 3: 1.146819e-39 < 0.1 TPM 3
## 4: 1.229960e-07 < 0.1 TPM 4
## 5: 3.585673e-16 < 0.1 TPM 5
## ---
## 519: 9.981573e-07 < 0.1 TPM 519
## 520: 1.812008e-04 < 0.1 TPM 520
## 521: 7.869110e-10 < 0.1 TPM 521
## 522: 6.365904e-15 < 0.1 TPM 522
## 523: 2.185946e-05 < 0.1 TPM 523
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: EDN1 143.61165 3.493408 0.3252202 10.406217 2.322691e-25
## 2: CCN2 3838.81312 3.311703 0.1858269 17.809254 5.990320e-71
## 3: GPR87 533.54706 3.282393 0.2390153 13.656904 1.836253e-42
## 4: ALPK2 44.93538 3.174662 0.4192791 5.808344 6.309395e-09
## 5: MUC12 19.79726 2.951746 0.4237826 6.511310 7.449834e-11
## ---
## 527: TRMT9B 145.42958 -1.727955 0.3060439 -5.633807 1.762745e-08
## 528: TBC1D3L 29.13473 -1.733263 0.3769001 -4.528622 5.936973e-06
## 529: CTXND2 171.45399 -1.825347 0.2641534 -6.887617 5.673462e-12
## 530: WFIKKN2 376.78945 -1.939575 0.2292372 -8.446372 3.004952e-17
## 531: MROH7-TTC4 20.24325 -2.216531 0.4251913 -4.950916 7.386496e-07
## padj FDR Input Rank
## 1: 8.487113e-23 < 0.1 Counts 1
## 2: 2.188863e-67 < 0.1 Counts 2
## 3: 1.118278e-39 < 0.1 Counts 3
## 4: 3.440975e-07 < 0.1 Counts 4
## 5: 5.917760e-09 < 0.1 Counts 5
## ---
## 527: 8.588092e-07 < 0.1 Counts 527
## 528: 1.583482e-04 < 0.1 Counts 528
## 529: 5.315597e-10 < 0.1 Counts 529
## 530: 4.773954e-15 < 0.1 Counts 530
## 531: 2.726289e-05 < 0.1 Counts 531
##
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: MROH7-TTC4 19.77094 -2.336210 0.4243946 -5.001629 5.684777e-07
## 2: WFIKKN2 371.00721 -1.937066 0.2298615 -8.412662 4.008097e-17
## 3: CTXND2 168.84610 -1.817134 0.2652655 -6.827901 8.616600e-12
## 4: TBC1D3L 28.55182 -1.742442 0.3813485 -4.499991 6.795649e-06
## 5: TRMT9B 145.07064 -1.716032 0.3049279 -5.610108 2.021999e-08
## ---
## 519: SERPINE1 8891.97515 2.850970 0.3263039 8.759388 1.963139e-18
## 520: MUC12 31.29951 3.050554 0.4135533 5.990845 2.087531e-09
## 521: GPR87 526.04952 3.283122 0.2391306 13.655049 1.883633e-42
## 522: CCN2 3782.15461 3.313824 0.1852169 17.879269 1.710657e-71
## 523: EDN1 141.58628 3.483585 0.3262539 10.347843 4.280208e-25
## padj FDR Input Rank
## 1: 2.185946e-05 < 0.1 TPM 1
## 2: 6.365904e-15 < 0.1 TPM 2
## 3: 7.869110e-10 < 0.1 TPM 3
## 4: 1.812008e-04 < 0.1 TPM 4
## 5: 9.981573e-07 < 0.1 TPM 5
## ---
## 519: 3.585673e-16 < 0.1 TPM 519
## 520: 1.229960e-07 < 0.1 TPM 520
## 521: 1.146819e-39 < 0.1 TPM 521
## 522: 6.249030e-68 < 0.1 TPM 522
## 523: 1.563560e-22 < 0.1 TPM 523
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: MROH7-TTC4 20.24325 -2.216531 0.4251913 -4.950916 7.386496e-07
## 2: WFIKKN2 376.78945 -1.939575 0.2292372 -8.446372 3.004952e-17
## 3: CTXND2 171.45399 -1.825347 0.2641534 -6.887617 5.673462e-12
## 4: TBC1D3L 29.13473 -1.733263 0.3769001 -4.528622 5.936973e-06
## 5: TRMT9B 145.42958 -1.727955 0.3060439 -5.633807 1.762745e-08
## ---
## 527: MUC12 19.79726 2.951746 0.4237826 6.511310 7.449834e-11
## 528: ALPK2 44.93538 3.174662 0.4192791 5.808344 6.309395e-09
## 529: GPR87 533.54706 3.282393 0.2390153 13.656904 1.836253e-42
## 530: CCN2 3838.81312 3.311703 0.1858269 17.809254 5.990320e-71
## 531: EDN1 143.61165 3.493408 0.3252202 10.406217 2.322691e-25
## padj FDR Input Rank
## 1: 2.726289e-05 < 0.1 Counts 1
## 2: 4.773954e-15 < 0.1 Counts 2
## 3: 5.315597e-10 < 0.1 Counts 3
## 4: 1.583482e-04 < 0.1 Counts 4
## 5: 8.588092e-07 < 0.1 Counts 5
## ---
## 527: 5.917760e-09 < 0.1 Counts 527
## 528: 3.440975e-07 < 0.1 Counts 528
## 529: 1.118278e-39 < 0.1 Counts 529
## 530: 2.188863e-67 < 0.1 Counts 530
## 531: 8.487113e-23 < 0.1 Counts 531
# Set a function rebuilding DE tables with gene ranks
combineTable_fn <- function(List){
# Select columns and join the data frames by gene
full_join(List[[TvC[1]]][,.(Gene, Input, Rank, baseMean)],
List[[TvC[2]]][,.(Gene, Input, Rank, baseMean)], by="Gene") %>%
# Add columns assining gene expression levels, rank differences, and mean ranks
mutate(logMeanExpression=log(baseMean.x+baseMean.y/2),
RankDiff=Rank.x-Rank.y,
MeanRank=(Rank.x+Rank.y)/2)
}
# Explore outputs of the function
head(combineTable_fn(FDRrankList))
## Gene Input.x Rank.x baseMean.x Input.y Rank.y baseMean.y
## 1: CCN2 TPM 1 3782.1546 Counts 1 3838.8131
## 2: F2R TPM 2 2460.4827 Counts 2 2498.1179
## 3: CDCP1 TPM 3 1500.9144 Counts 3 1523.3025
## 4: STK38L TPM 4 11846.6666 Counts 4 12016.7778
## 5: TM4SF18 TPM 5 1895.6825 Counts 5 1919.1263
## 6: GPR87 TPM 6 526.0495 Counts 6 533.5471
## logMeanExpression RankDiff MeanRank
## 1: 8.648495 0 1
## 2: 8.218664 0 2
## 3: 7.724255 0 3
## 4: 9.790042 0 4
## 5: 7.956913 0 5
## 6: 6.675600 0 6
dim(combineTable_fn(FDRrankList))
## [1] 535 10
tail(combineTable_fn(FDRrankList))
## Gene Input.x Rank.x baseMean.x Input.y Rank.y baseMean.y logMeanExpression
## 1: UTP3 <NA> NA NA Counts 524 308.86592 NA
## 2: IGIP <NA> NA NA Counts 525 90.33876 NA
## 3: CUTALP <NA> NA NA Counts 527 80.56840 NA
## 4: CREG1 <NA> NA NA Counts 529 273.30711 NA
## 5: P2RY13 <NA> NA NA Counts 530 74.69163 NA
## 6: RBM43 <NA> NA NA Counts 531 33.32318 NA
## RankDiff MeanRank
## 1: NA NA
## 2: NA NA
## 3: NA NA
## 4: NA NA
## 5: NA NA
## 6: NA NA
# Set a function plotting gene ranks between TPM- (x-axis) and Counts-Inputs (y-axis)
compareRanks_fn <- function(df, rankedby) {
ggplot(df,
aes(x=Rank.x,
y=Rank.y,
color=logMeanExpression)) +
geom_point(alpha=0.5) +
theme_bw() +
xlab("Rank with TPM") +
ylab("Rank with Counts") +
geom_abline(slope=1, color="black", size=0.5) +
ggtitle(paste(rankedby, "Ranking with TPM vs Counts Inputs")) +
scale_color_gradient(low="blue", high="red")
}
# Set a function plotting the rank difference over the gene expression level
RankdiffOverMean_fn <- function(df, rankedby) {
ggplot(df,
aes(x=logMeanExpression,
y=RankDiff,
color=MeanRank)) +
geom_point(alpha=0.5) +
theme_bw() +
ylab("Rank Difference (TPM - Counts)") +
ggtitle(paste("Rank Difference Inputs (TPM - Counts) in", rankedby)) +
geom_hline(yintercept=0, color="black", size=0.5) +
scale_color_gradient(low="blue", high="red")
}
# Ranked by FDR
compareRanks_fn(combineTable_fn(FDRrankList),
"FDR")
# Ranked by absolute fold change
compareRanks_fn(combineTable_fn(lfcList),
"Absolute Log2FoldChange")
# Ranked by magnitude of positive fold change
compareRanks_fn(combineTable_fn(UPlfcList),
"Log2FoldChange (Increased)")
# Ranked by magnitude of negative fold change
compareRanks_fn(combineTable_fn(DOWNlfcList),
"Log2FoldChange (Decreased)")
# Ranked by FDR
RankdiffOverMean_fn(combineTable_fn(FDRrankList),
"FDR")
# Ranked by absolute fold change
RankdiffOverMean_fn(combineTable_fn(lfcList),
"Absolute Log2FoldChange")
# Ranked by magnitude of positive fold change
RankdiffOverMean_fn(combineTable_fn(UPlfcList),
"Log2FoldChange (Increased)")
# Ranked by magnitude of negative fold change
RankdiffOverMean_fn(combineTable_fn(DOWNlfcList),
"Log2FoldChange (Decreased)")
# Clean data to generate an upset plot
res <- res %>%
# Filter genes with valid padj
filter(!is.na(padj)) %>%
# Detect genes which are up/down/unchanged change patterns in either TPM and Counts inputs
mutate(Up=ifelse(FDR == FDRv[1] & log2FoldChange > 0, Gene, ""), # What are upregulated genes?
Down=ifelse(FDR == FDRv[1] & log2FoldChange < 0, Gene, ""), # What are downregulated genes?
Unchanged=ifelse(FDR == FDRv[2], Gene, ""), # What are unchanged genes?
TPM_Input=ifelse(Input == "TPM", Gene, ""), # What are the genes from TPM input?
Counts_Input=ifelse(Input == "Counts", Gene, "")) # What are the genes from Counts input?
# Create a list storing groups of interest
upsetInput <- list(Up=res$Up,
Down=res$Down,
Unchanged=res$Unchanged,
TPM_Input=res$TPM,
Counts_Input=res$Counts)
# Create the upset plot
upset(fromList(upsetInput),
sets.x.label="Number of Genes per Group",
order.by="freq")
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS
##
## Matrix products: default
## BLAS/LAPACK: /home/mira/miniconda3/envs/r/lib/libopenblasp-r0.3.10.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] AnnotationDbi_1.52.0 UpSetR_1.4.0
## [3] ggrepel_0.8.2 ggplotify_0.0.5
## [5] gridExtra_2.3 pheatmap_1.0.12
## [7] DESeq2_1.30.0 SummarizedExperiment_1.20.0
## [9] Biobase_2.50.0 MatrixGenerics_1.2.0
## [11] matrixStats_0.57.0 GenomicRanges_1.42.0
## [13] GenomeInfoDb_1.26.2 IRanges_2.24.0
## [15] S4Vectors_0.28.1 tximport_1.18.0
## [17] forcats_0.5.0 stringr_1.4.0
## [19] dplyr_1.0.2 purrr_0.3.4
## [21] readr_1.4.0 tidyr_1.1.2
## [23] tibble_3.0.4 ggplot2_3.3.2
## [25] tidyverse_1.3.0 AnnotationHub_2.22.0
## [27] BiocFileCache_1.14.0 dbplyr_2.0.0
## [29] BiocGenerics_0.36.0 rmarkdown_2.5
## [31] data.table_1.13.4
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-0 ellipsis_0.3.1
## [3] XVector_0.30.0 fs_1.5.0
## [5] rstudioapi_0.13 farver_2.0.3
## [7] bit64_4.0.5 interactiveDisplayBase_1.28.0
## [9] fansi_0.4.1 lubridate_1.7.9.2
## [11] xml2_1.3.2 splines_4.0.3
## [13] geneplotter_1.68.0 knitr_1.30
## [15] jsonlite_1.7.2 broom_0.7.2
## [17] annotate_1.68.0 shiny_1.5.0
## [19] BiocManager_1.30.10 compiler_4.0.3
## [21] httr_1.4.2 rvcheck_0.1.8
## [23] backports_1.2.1 assertthat_0.2.1
## [25] Matrix_1.2-18 fastmap_1.0.1
## [27] cli_2.2.0 later_1.1.0.1
## [29] htmltools_0.5.0 tools_4.0.3
## [31] gtable_0.3.0 glue_1.4.2
## [33] GenomeInfoDbData_1.2.4 rappdirs_0.3.1
## [35] Rcpp_1.0.5 cellranger_1.1.0
## [37] vctrs_0.3.5 xfun_0.19
## [39] ps_1.5.0 rvest_0.3.6
## [41] mime_0.9 lifecycle_0.2.0
## [43] XML_3.99-0.5 zlibbioc_1.36.0
## [45] scales_1.1.1 hms_0.5.3
## [47] promises_1.1.1 RColorBrewer_1.1-2
## [49] yaml_2.2.1 curl_4.3
## [51] memoise_1.1.0 stringi_1.5.3
## [53] RSQLite_2.2.1 BiocVersion_3.12.0
## [55] genefilter_1.72.0 BiocParallel_1.24.1
## [57] rlang_0.4.9 pkgconfig_2.0.3
## [59] bitops_1.0-6 evaluate_0.14
## [61] lattice_0.20-41 labeling_0.4.2
## [63] bit_4.0.4 tidyselect_1.1.0
## [65] plyr_1.8.6 magrittr_2.0.1
## [67] R6_2.5.0 generics_0.1.0
## [69] DelayedArray_0.16.0 DBI_1.1.0
## [71] pillar_1.4.7 haven_2.3.1
## [73] withr_2.3.0 survival_3.2-7
## [75] RCurl_1.98-1.2 modelr_0.1.8
## [77] crayon_1.3.4 locfit_1.5-9.4
## [79] grid_4.0.3 readxl_1.3.1
## [81] blob_1.2.1 reprex_0.3.0
## [83] digest_0.6.27 xtable_1.8-4
## [85] httpuv_1.5.4 gridGraphics_0.5-0
## [87] munsell_0.5.0